package com.numericalmethod.algoquant.model.ralph2009.draft;

import java.util.Random;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math3.linear.RealMatrix;

/**
 * 
 * Implementation of Single Shock dynamics in section 4 of the paper.
 * 
 * 
 */
public class SingleShockPriceDynamics extends PriceDynamics {

	// parameters to control the shock
	private int shockPeriodStartStep;
	private int shockPeriodEndStep;
	private boolean shockWithGaussian;
	
	public SingleShockPriceDynamics(double mu0, double mu1, double sigmaS1, double alpha, double sigmaX1, double sigmaX2, double phi, 
			int shockPeriodStartStep, int shockPeriodEndStep, boolean shockWithGaussian) {
		super(mu0, mu1, sigmaS1, alpha, sigmaX1, sigmaX2, phi);

		this.shockPeriodStartStep=shockPeriodStartStep;
		this.shockPeriodEndStep=shockPeriodEndStep;
		this.shockWithGaussian=shockWithGaussian;
	}
	
	@Override
	public void simulatePath(){
		//assume stock price follows log normal
		Random rand1=new Random(2965179875022589598L);
		
		Random rand2=new Random(17399225432l);
		
		this.z1Matrix = new double[pathCount][timePeriod];
		this.z2Matrix = new double[pathCount][timePeriod];

		
		for (int i = 0; i < pathCount; i++) {
			for (int j = 0; j < timePeriod; j++) {
				double r1 = rand1.nextGaussian();
				double r2 = rand2.nextGaussian();
				// only shock once (According to Section 4 of the paper)
				if (j >= shockPeriodStartStep && j <= shockPeriodEndStep) {

					if (shockWithGaussian) {
						// full randomness
						z1Matrix[i][j] = r1;
						z2Matrix[i][j] = r2;
					}
					else {
						// 1.0 shock but allow a very little randomness otherwise
						// regression will break down.
						z1Matrix[i][j] = 1.0 + r1 * Math.sqrt(0.000001);
						z2Matrix[i][j] = 1.0 + r2 * Math.sqrt(0.000001);
					}

				} else {
					// no shock if it is not within the "shocking time"
					/*
					 * CANNOT SIMPLY PUT THEM TO ZEROs BECAUSE THAT WOULD FAIL
					 * THE REGRESSION at regATPlus1Vec =
					 * LSRegression.estimateYValues(excessReturnVector,
					 * divYieldVector, aTPlus1Vec); z1Matrix[i][j]=0.0;
					 * z2Matrix[i][j]=0.0;
					 */
					// Suppress the randomness to very small variance, such that
					// they almost zero movements
					// but still allow the regression to pass without error.
					z1Matrix[i][j] = r1 * Math.sqrt(0.000001);
					z2Matrix[i][j] = r2 * Math.sqrt(0.000001);

					// z1Matrix[i][j]=r1;
					// z2Matrix[i][j]=r2;

				}
			}
		}
		
		runDividendSimulation();
		runMomentumSimulation();
		runStockSimulation();
	
	}	

	//The S(t) simulation, i.e. Eq. 13
	// dSt = St[ ((mu0 + mu1*Xt)(1-phi) + phi*Mt)*dt + sigmaS1*dZt ] 
	@Override
	protected void runStockSimulation(){
		if(stockPriceMatrix==null){
			runMomentumSimulation();
			final RealMatrix localMomentumMatrix = this.momentumMatrix;
			final RealMatrix localDivYieldMatrix = this.divYieldMatrix;
			stockPriceMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			excessReturnMatrix = new Array2DRowRealMatrix(pathCount, timePeriod-1);
			stockPriceMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return s0;
					}else{
						double z1 = z1Matrix[row][column];
						double lastMt = localMomentumMatrix.getEntry(row, column-1);
						double lastXt = localDivYieldMatrix.getEntry(row, column-1);
						double lastSt = stockPriceMatrix.getEntry(row, column-1);
						
						
						// dS/S = (...)Xt + (...)Mt + AR_coef * dS(t-1)/S(t-1) + AR_coef * lastMt + AR_coef * lastXt + innovation
						double dSt = lastSt*( ((mu0 + mu1*lastXt)*(1-phi) + phi*lastMt)*dt + sigmaS1*Math.sqrt(dt)*z1 );
						
						excessReturnMatrix.setEntry(row, column-1, (dSt/lastSt) - (riskFree*dt));
						return lastSt + dSt;
					}
				}
			});
		}
	}
	
	//The M(t) simulation, i.e. Eq. 14 
	// dMt = (1 - phi)(mu0 + mu1*Xt - Mt*) *dt + sigmaS1 * dZt
	@Override
	protected void runMomentumSimulation(){
		if(momentumMatrix == null){
			momentumMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			runDividendSimulation();
			final RealMatrix localDivYieldMatrix = this.divYieldMatrix;
			momentumMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return m0;
					}else{
						double z1 = z1Matrix[row][column];
						
						
						double lastMt = momentumMatrix.getEntry(row, column-1);
						double lastXt = localDivYieldMatrix.getEntry(row, column-1);
						double dMt = (1-phi)*(mu0+ (mu1*lastXt) - lastMt) *dt
								+ (sigmaS1 * Math.sqrt(dt) * z1);
						return lastMt + dMt;
					}
				}
			});	
		}
	}

	
	//The X(t) simulation, i.e. Eq. 15 (Xt follows OU process)
	// dXt = -alpha * Xt * dt + sigmaX' * dZt and 
	// mut = mu0 + mu1 * Xt
	// or, 
	// Xt+1 = Xt + dXt
	@Override
	protected void runDividendSimulation(){
		if(divYieldMatrix == null){
			divYieldMatrix = new Array2DRowRealMatrix(pathCount, timePeriod);
			divYieldMatrix.walkInColumnOrder(new DefaultRealMatrixChangingVisitor(){
				public double visit(int row, int column, double value){
					if(column == 0){
						return x0;
					}else{
						double z1 = z1Matrix[row][column];
						double z2 = z2Matrix[row][column];
						

						
						double lastXt = divYieldMatrix.getEntry(row, column-1);
						
						double dXt = (-alpha * dt * lastXt)

								+ (sigmaX1*Math.sqrt(dt)*z1) +(sigmaX2*Math.sqrt(dt)*z2);
						
						return lastXt + dXt;
					}
				}
			});		
		}
	}	
	

}
